pool_pred_summary <- function(list_pred_summaries) {
  ## Sub-functions
  make_row_vector <- function (x) {
    if (is.vector(x) & !is.matrix(x)) {
      return(t(as.matrix(x)))
    } else {
      return(x)
    }
  }
  
  pool_imputations <- function(list_pred_summaries) {
    est <- sapply(list_pred_summaries, function(x)
      x$Prediction) %>%
      make_row_vector()
    se <- sapply(list_pred_summaries, function(x)
      x$SE) %>%
      make_row_vector()
    n <- sapply(list_pred_summaries, function(x)
      x$n) %>%
      make_row_vector()
    sd <- sapply(list_pred_summaries, function(x)
      x$sd) %>%
      make_row_vector()
    
    ## Point estimates
    est_pooled <-  apply(est, 1, mean)
    
    ## Standard errors
    pool_ses <- function(est, se) {
      m <- ncol(se)
      V_W <- apply(se, 1, function(x)
        mean(x ^ 2))
      V_B <- (1 - m) * apply(est, 1, function (x)
        sum(x - mean(x)))
      V_T <- V_W + V_B + V_B / m
      SE_pooled <- sqrt(V_T)
      return(SE_pooled)
    }
    ses_pooled <- pool_ses(est, se)
    
    ## n, sd
    n_pooled <- apply(n, 1, mean)
    sd_pooled <- apply(sd, 1, mean)
    
    ## Value
    return(data.frame(Prediction = est_pooled,
                      SE = ses_pooled,
                      n = n_pooled,
                      sd = sd_pooled))
  }
  
  ## Pool estimates
  pooled_imputations <- pool_imputations(list_pred_summaries)
  
  ## Initialize ourput object
  pooled_summary <- list_pred_summaries[[1]]
  
  ## Drop unneccesary columns
  pooled_summary$Prediction <-
    pooled_summary$SE <-
    pooled_summary$z <-
    pooled_summary$p <-
    pooled_summary$lower <-
    pooled_summary$upper <- 
    NULL
  
  ## Pool point estimates and SEs
  pooled_summary$Prediction <- pooled_imputations$Prediction
  pooled_summary$SE <- pooled_imputations$SE
  pooled_summary$n <- pooled_imputations$n
  pooled_summary$sd <- pooled_imputations$sd
  
  ## Add upper and lower bounds
  pooled_summary$lower <-
    pooled_summary$Prediction + qnorm(.025) * pooled_summary$SE
  pooled_summary$upper <-
    pooled_summary$Prediction + qnorm(.975) * pooled_summary$SE
  
  
  ## Return
  return(pooled_summary)
}